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The study on the coUective dynamics of synchronization among genetic osciUators is essential for 
the understanding of the rhythmic phenomena of hving organisms at both molecular and cellular 
levels. Genetic oscillators are biochemical networks, which can generally be modelled as nonlinear 
dynamic systems. We show in this paper that many genetic oscillators can be transformed into Lur'e 
form by exploiting the special structure of biological systems. By using control theory approach, 
we provide a theoretical method for analyzing the synchronization of coupled nonidentical genetic 
oscillators. Sufficient conditions for the synchronization as well as the estimation of the bound of the 
synchronization error are also obtained. To demonstrate the effectiveness of our theoretical results, 
a population of genetic oscillators based on the Goodwin model are adopted as numerical examples. 

PACS numbers: 87.16.Yc, 05.45.Xt, 89.75.Hc 



I. INTRODUCTION 

Elucidating cooperative behavior of synchronization of coupled genetic oscillators has important biological implica- 
tions and potential engineering applications from both theoretical and experimental viewpoints, and it is also essential 
for the understanding of the rhythmic phenomena of living organisms at both molecular and cellular levels. So far 
many researchers have studied the synchronization in genetic networks from the aspects of experiment, numerical 
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simulation and theoretical analysis. For instance, in ^, the authors experimentally investigated the synchronization 
of cellular clock in the suprachiasmatic nucleus (SCN); in y and y, the collective dynamics of synchronization are 
theoretically studied in synthetic biological networks of identical genetic oscillators; and in (^], the mechanism of syn- 
chronization in a population of identical hysteretic genetic oscillators is analyzed. Biologically, the genetic oscillators, 
even in the same species, are usually nonidentical possibly due to asymmetrical nutrition conditions and fluctuated 
environments, and the nonidentical property can be modelled as parametric mismatches among oscillators. For ex- 
ample, in the suprachiasmatic nucleus (SCN), the periods of the circadian oscillators are not exactly the same, and 
it has been observed that isolated individual neurons are able to produce circadian oscillations with period ranging 
from 20 to 28 hours 0, 3l ■ y 1 the synchronization for nonidentical genetic oscillators is examined numerically. 
Although many mathematical models have been developed to study the cooperative behaviors of cellular oscillation, 
there is no general theoretical method in analyzing the dynamics of the coupled genetic oscillators due to their inherent 
nonlinearity. 

Genetic networks are biochemically dynamical systems, in which the nodes indicate the biochemicals, and the 
couplings represent the biochemical interactions llo| . Mathematically many genetic oscillators can be expressed in 



al svst( 



the form of multiple additive terms, each of which particularly is of linear, Michaelis-Menten or Hill forms, such as 
the well-known Goodwin model repressilator 13,j, toggle switch 14], and the circadian oscillators [l5|. From 
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the synthetic biology viewpoint, genetic oscillators with only linear, Michaelis-Menten and Hill form terms can also 
be implemented easily. In this paper we explore such special structure of gene networks to show that these genetic 
oscillators can be transformed into Lur'e form and can be further analyzed by using Lur'e system method in control 
theory 

The aim of this paper is to provide a general theoretical method for analyzing the synchronization of coupled 
nonidentical genetic oscillators with the above-mentioned structure. In studying the synchronization of genetic oscil- 
lators (and other nonlinear systems), a general idea is to study the stability of the error equations among oscillators. 
However, there are two main difficulties for such method: (1) the difference of the oscillator dynamics usually cannot 
be written into a function of the state error; (2) due to the nonlinearity, there is no general efficient analysis method 
for the stability of the error system. We show in this paper that for genetic oscillators with the above-mentioned 
structure, we can overcome both of the above difficulties. Since coupled nonidentical oscillators usually cannot achieve 
complete synchronization, a synchronous error is required to evaluate the quality of the synchronization. In this paper, 
we present a theoretical result, which can not only guarantee the synchronization, but also estimates the bound of 



the synchronization error. Besides, the obtained conditions can be represented in terms of Unear matrix inequahties 



(LMIs] 



which are very easy to be verified. Recenth^it was found that many biological networks are complex 



networks with small-world and scale-free properties 



cently, r 



Our method is also applicable to genetic oscillator net- 
works with complex topology, directed and weighted couplings. To demonstrate the effectiveness of the theoretical 
results, we present two simulation examples of coupled Goodwin oscillators with linear and Michaelis-Menten cou- 
plings, respectively. Finally, several remarks on the extensions of the proposed method are discussed. Notation used 
in this paper as well as the detailed theoretical analysis are given in the Appendix. 

II. METHODS AND RESULTS 

Mathematical modelling provides a powerful tool for studying gene regulation processes in living organisms. Basi- 
cally, there are two types of genetic network models, i.e., the Boolean model (or discrete model) and the differential 



equation model (or continuous model) 
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In Boolean models, the activity of each gene is expressed in one 
of two states, ON or OFF, and the state of a gene is determined by a Boolean function of the states of other related 
genes. In the differential equation models, the variables describe the concentrations of gene products, such as mRNAs 
and proteins, as continuous values, which are more accurate and can provide detailed understanding of the dynamical 
behaviors of the gene regulation systems. In this paper, by adopting the differential equation models, we consider the 
following form of a general genetic oscillator: 

I 

1=1 

where y{t) £ i?" represents the concentrations of proteins, RNAs and chemical complexes, A and Bi are matrices in 
-R"^", fi{y) — [fiiiuiit)), • • • , fin{yn{t))]'^ with fij{yj{t)) as a monotonic increasing or decreasing regulatory function, 
which usually is of the Michaelis-Menten or Hill form. Undoubtedly, many well-known genetic system models can be 
represented in this form, such as the Goodwin model 



, and the 
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[111, [131, the repressilator ISj, the toggle switch |14j 
circadian oscillators ^isj . In synthetic biology, genetic oscillators of this form can be implemented experimentally [3 ■ 
To make our method more understandable and to avoid unnecessarily complicated notation, we consider the following 
simplified model, in which there are only one increasing and one decreasing nonlinear terms in each equation of the 
genetic oscillator. 

m = Ay{t) + B^h{y{t)) -f 52/2(?/(t)), (2) 



where Ay{t) includes the degradation terms and all the other linear terms in the genetic oscillator, fi{y{t)) = 
[fii{yi{t)), ■ ■ ■ , fin{yn{t))]'^ with fij{yj{t)) as a monotonic increasing function of the Hill form 

and f2{y{t)) = [/2i(j/i(i)), • • • , f2n{yn{t))]'^ with f2j{vj{t)) as a monotonic decreasing function of the following form 

In the above equations, both Hij and iJ2j are the Hill coefficients. To avoid confusion, we let the jth column of B12 
be zeros if fij,2j = 0. Since 

by letting /(.) = /i(.), we can rewrite (2) as follows: 

y{t) = Ay{t) + B,f{y{t)) - B2g{y{t)) + B2. (3) 

Obviously, fi and gi satisfy the following sector conditions 

n < h{a)-h{b) < -L . 

n < g^{a■)-gi{b) ^ , 
" - a-f, - '^^Sj, 



Va, 6 € i? (a 7^ 6); i = 1, • ■ • , n, 



or equivalently. 



(/.(a) - /»(^^))[(/^(a) - m) - ku{a - b)] < 0, 

(5. (a) - ff,;(&))[(g.(a) - - fc2.(a - &)] < 0, (4) 

Va, 6 G i? (a 7^ 6); i = 1, • • • , n, 
It follows from the mean value theorem that for differentiable fi and gi, the above sector conditions correspond to 

0<f (a)<fcH, 

< ^(a) < fc2., (5) 
Va e i?; i = 1, ■ • ■ , n. 

Recall that a Lur'e system is a linear dynamic system, feedback interconnected to a static nonlinearity /(.) that 

n 

satisfies a sector condition Ifjj. Hence, the genetic oscillator (3) can be seen as a Lur'e system, which can be 
investigated by using the fruitful Lur'e system method in control theory. In the following, we first consider coupled 
identical genetic oscillators, and then extend the result to the nonidentical case. 
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We first analyze the following N linearly coupled genetic oscillators, in which each genetic oscillator is identical. 

N 

x,{t) = Ax,{t) + Bif{x,{t)) - B2g{x,{t)) + B2+J2 G^iDxj{t), z = 1, • • • , iV (6) 

where Xi{t) G i?" is the state vector of the ith genetic oscillator (corresponds to y{t) in Eq. (3)), D G defines 
the coupling between two genetic oscillators. G = {Gij)NxN is the coupling matrix of the network, in which Gij is 
defined as follows: if there is a link from oscillator j to oscillator i {j ^ i), then Gij equals to a positive constant 
denoting the coupling strength of this link; otherwise, Gij — 0; Gu = Matrix G defines the coupling 

topology, direction, and the coupling strength of the network. 

Since in biological networks, the genetic oscillators are usually nonidentical, there are parametric mismatches among 
oscillators. Next, we consider the following network of N coupled nonidentical genetic oscillators: 

x^{t) ^{A + AA^{t))x,{t) + {Bi + ABu{t))f{x^{t)) - {B2 + AB2^{t))g{x^{t)) + {B2 + AB2{t)) 

(7) 

+ j:J=iG.jDx,{t),i^l,--- ,N 
where AAi, ABu, AB2i are the mismatch matrices, which can be time- varying. We assume that the mismatch matrices 
AAiit) , ABiiit) , AB2i{t) can be estimated by the following bounds, which are also reasonable for general biological 
systems. 

\\AA,lt)\\ < ai, \\ABu(t)\\ < «2, ||AB2.(t)|| < a^yi. 

We also assume that 

\\x^m < ^1, \\f{x,m\ < 62, \\g{x,m\ < ^3,vz. 

Since in genetic oscillators, Xi{t) usually denotes the concentrations of mRNA, protein, neurotransmitter, etc., which 
are of limited values, and /(.) and g{.) are usually monotonic functions with saturated values, the above assumptions 
are also reasonable. The other parameters are defined as the same as those in the identical case. 

For the above two network, i.e. Eq. (6) and (7), we mainly use Wu's method to analyze the synchronization 
which can separate the effects of the coupling and the individual genetic oscillator dynamics. Based on Lyapunov 
method, we can obtain the following sufficient conditions for the synchronization of coupled identical genetic oscillators 
(Eq. (6)) and coupled nonidentical genetic oscillators (Eq. (7)) in Theorems 1 and 2, respectively. In the following 
theorems and hereafter, Ki — diag(fcii, • • • , fci„), K2 — diag(fc2i, • • • , fc2„), and matrix U £ ji^xN defined as an 
irreducible matrices with zero row sums, whose off-diagonal elements are all non-positive. \min{P) and Xmax{P) 
represent the minimal and maximal eigenvalues of the matrix P respectively, and (E5 indicates the Kronecker product, 
which is defined in Appendix. 



(8) 



Theorem 1: If there are matrices P > 0, Ai = diag(Aii, • • • , Ai„) > 0, A2 = diag(A2i, • • • , A2„) > 0, Q e ii"^", 
and U G R^^^ , such that the foUowing matrix inequaUties hold 

PA + A^P-Q-Q'^ PBi+KiM -PB2 + K2A2 
Mi= BfP + KiAi -2Ai <0, 

-SjP + JsTsAa -2A2 

{UG ® PD + U ®QY + {UG ®PD + U ®Q)<Q, 
then the network (6) is asymptoticaUy synchronous. 

Theorem 2: If there are matrices P > 0, Ai = diag(Aii, • • • , Ai„) > 0, A2 = diag(A2i, • • • , A2„) > 0, Q G ii"^", 
U G R^^^ , and a positive real constant 7 such that the following matrix inequalities hold 

P^ + ^^P-g-Q^ + 7/ PBi+KiAi -PB2+K2A2 
M2= BfP + KiAi -2Ai <0, 

-BjP + K2A2 -2A2 

{UG (g) PD + U (g) + {UG (g) PL» + [/ (g) Q) < 0, 
then the network (7) is asymptotically synchronous with error bound 



(9) 



Y^{-u,,)\\^.{t) - x,{t)r < ^Y^E(-u^J) 



where 



raax 



The proof of the above theorems are somewhat technical, and we defer the details to the Appendix. In Theorem 2, 
we not only give a sufficient condition for the synchronization, but also provide an estimation of the synchronization 
error bound. In this error bound estimation, we can select the form of the matrix U beforehand to obtain different 
error combinations. Specifically, if we choose the following form of U 

N-1 -1 ... -1 
-1 1 

U = 

-1 

we can obtain the following synchronous error estimation: 



J2\\^j{t)-xr{t)f<{N-l) 

j=2 



max 
7 Xmin 



It should be note that, since in each step, we used conservative estimations of the bounds, the estimated error bound 
may be much larger than the actual error. In other words, if accurate information on the parametric mismatches is 
known, we can have a better error estimation. 

The first matrix inequality in (8) (or (9)) is an LMI, which can be easily verified by using convex optimization 
techniques, e.g., the interior point method [17 1, and by software packages, e.g., the MATLAB LMI Toolbox. If we 
choose the matrix U beforehand, the second matrix inequality in (8) (or (9)) is also an LMI. Furthermore, for two 
special cases, we have the following results 2^: 

(1) When G is symmetric: Letting U ~ —G, the second matrix inequality in (8) (or (9)) can be rewritten as 

-G^ ® {PD + D'^P) - G ® (Q + Q^) < 0, 

which is equivalent to 

ai{PD + D^P) + {Q + Q^) < 0, for all nonzero eigenvalues Ui of G. 

(2) When D is symmetric and commutable with G: Letting U = — (G + G^)/2, the second matrix inequality in (8) 
(or (9)) can be rewritten as 

-i(G + G^f (^PD-'^{G + G'^)®{Q + Q^) < 0, 

which is equivalent to 

2(JiPD + (Q + Q'^) < 0, for all nonzero real part cr,; of the eigenvalues of G. 

Thus the second matrix inequality in (8) (or (9)) is also a lower dimensional LMI, which together with the first 
LMI can be verified easily. 

Next, we use numerical examples to show the effectiveness of the theoretical results. 

III. NUMERICAL EXAMPLES 

To demonstrate the effectiveness of our theoretical methods, we study a population of coupled SCN neuron model 
oscillators. The single cell or genetic oscillator is described by the classical Goodwin model [ll|. In this model, a 
clock gene mRNA (A) produces a clock protein (Y), which activates a transcriptional inhibitor (Z). Z inhibits the 
transcription of the clock gene, thus forming a negative feedback loop. In this paper, we assume that the light L — 0. 
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The oscillators coupled through the release and receiving of neurotransmitter among neurons. Similar to j^, the 
evolution equations for a network of N coupled nonidentical oscillators are given below: 

Xi =^ Vii - Vi2Xi + KF, 

Yi = vaXi - ViiYi, 

(10) 

Vi = VitXi - Vi^Vi, 

where Vii^Vi2TVa,Viii,Vi^,Vi(,,Vi'j,Vi^ are positive constants, and > is the coupling strength. The variables 
Xi,Yi^ Zi describe the dynamics of the oscillator in the zth neuron, and Vi describes the evolution of the neurotrans- 
mitter in the ith neuron. The release of the neurotransmitter is supposed to be fast with respect to the 24-h timescale 
of the oscillators and becomes homogeneous to establish an average neurotransmitter level, or a mean field F 

1 ^ 
N ^ 

Clearly, the individual Goodwin model is of the form (3), in which / = 0, Bi = 0, 5 = [0, 0, + ), 0]^, 

is a 4 X 4 matrix with all zero entries except for -62(1, 3) = ua, and all the other terms are in the linear form. By plus 
and minus KVi in the first equation of (10), we can get the coupling term ~ ^i)- 

The purpose of this section is to demonstrate the effectiveness of the theoretical method, instead of mimicking 
the real SCNs. We consider a small size of network with = 10 coupled oscillators, although there are ^ lO'' 
neurons in the SCNs. The concentrations are expressed in riM, and the standard parameters are set as m = 12, Vi\ — 
InM /h,Vi^ — Vi5 = l/h,Vii — Q.2/h,Vi2 = 0.3/h,Vi4 — Q.22/h,ViQ = Q.15/h,Vis — 2/h for all i, so that the period of 
the oscillator is approximately 24h, and the coupling strength is K ^ 1. It is known that the period of the Goodwin 
model is sensitive to the parameters Vi2,Vi4,Vie jl2]. The mismatches of Vi2,Vi4,Vie are randomly distributed in 
±10% around the above values of Vi2,Vi4,ViG. In Fig. 1, when starting from the same initial values, we show the 
oscillation dynamics of the niRNA concentrations of the 10 uncoupled oscillators, which indicate that the periods 
of the oscillators are quite different. Submitting the above parameters to the corresponding matrices in the matrix 
inequalities (9) of Theorem 2, letting U = — G and using MATLAB LMI Toolbox, we can easily find feasible solutions 
for (9), which indicate that the above all-to-all coupled network can achieve synchronization although it is not a 
complete synchronous state. In Fig. 2 (a), when starting from different initial values, we plot the time evolution of 
the mRNA concentrations (Xi) of all the oscillators. Fig. 2 (b) shows the synchronization error 

N 

J = Y.KX, - Xi)2 + (Y, - yi)2 + (z, - Zif + (V, - Vif]. 

i=2 
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which is gradually reduced with time evolution. 
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FIG. 1: Oscillation dynamics of the mRNA concentrations of the uncoupled oscillators with the same initial conditions. 

Since there is a maximal activity of fully active promoters, in 0|, the authors considered a Michaelis-Menten form 
of coupling term, that is, in (10) replacing KF by K^^. Our theoretical results in Theorems 1 and 2 are also 
applicable for this case. Specifically, we may have different methods to treat the coupling term. One of the simplest 
ways is explained as follows: Since the coupling term is the same for all oscillators and U is zero row sums, it 

is easy to show that the product oi U (E) P and the vector of the coupling terms (containing the coupling terms of all 
the oscillators) is zero. By using the same analysis as that of Theorem 2, it is easy to show that the synchronization 
condition is only the first LMI in Theorem 2 with Q — 0. For the Michaelis-Menten coupling, Figs. 3 (a) and 3 
(b) show the time evolution of the mRNA concentrations and the synchronization error, respectively. Figs. 2 and 3 
indicate that the coupled oscillators are indeed synchronized with small error bounds, which confirms the theoretical 
results. 



IV. CONCLUSION AND OUTLOOK 



In this paper, we presented a theoretical method for analyzing the synchronization of coupled nonidentical genetic 
oscillators based on control theory approach. The purpose of this paper is not to provide a general theory for all 
genetic oscillator networks, but provide an efficient method for genetic oscillator networks that can be expressed 
in the form of (1). In addition, the sufficient conditions for the synchronization were also derived based on LMI 
formalism, which can be easily verified numerically. Although the method is proposed for genetic oscillator networks, 
it is also applicable to other biochemical and neuronal networks of the form. To make the theoretical method more 
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FIG. 2: Time evolution of Goodwin models with linear coupling. (a)Tlie time evolution of the mRNA concentrations of the 
oscillators; (b) The time evolution of the synchronization error. 

understandable and to avoid unnecessarily complicated notation, we discussed only on some simplified forms of the 
genetic oscillators, but more general cases and extensions regarding this topic can be studied in a similar way, for 
example: 

1. The theoretical results can be easily extended to the general form of genetic oscillators with more than 2 
nonlinear terms in each equation as shown in (1). 

2. The genetic oscillator model (3) can be extended to a more general case such that fi and g,;, the components of 
f{y{t)) and g{y{t)), are functions oiy{t), instead of yi{t). For this case, we only require that 



< 



Va,6 e i?" (a ^b)]i = 1, 
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FIG. 3: Time evolution of Goodwin models with Michaelis-Menten coupling. (a)Tlie time evolution of the mRNA concentrations 
of the oscillators; (b) The time evolution of the synchronization error. 

where Cii,C2i G i?" are arbitrary nonzero real vectors. Moreover, fi and gi can be of more complex forms 
(non-Hill form) and non-differentiable provided that they satisfy sector conditions. 

3. Genetic oscillators who have a few terms that are not in linear, Michaelis-Menten and Hill forms, can also be 
analyzed similarly by using our method. For example, in the mammalian circadian clock model [2y|, there are a 
few product terms (of the form xy) besides the three kinds of terms. We can treat these terms as follows: Since 

xivi - X2y2 = xi(yi - y2) + y2[xi - X2), and < max{xi), < max(y2), 

yi - 2/2 Xi- X2 

then we can treat these terms similarly as those that satisfies the sector condition. Although such manipulation 
introduces additional conservation, the conservation is assumed to be limited since most terms in the model are 
of the three kind of terms. 
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|3l|, 



In addition to the nonlinear and coupling properties, genetic networks are intrinsically noisy |27 
and with significant time delays 

1313 

Future works regarding this topic also include the extension of our method 
to the case with noise perturbations and time delays. 
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GLOSSARY 

Gene regulatory network or genetic network: A gene regulatory network is a collection of DNA segments in a cell 
which interact with each other and with other substances in the cell, thereby governing the rates at which genes in 
the network are transcribed into mRNA. 

Circadian rhythm: Circadian rhythm is the name given to the roughly 24 hour cycles shown by physiological 
processes in plants, animals, fungi and cyanobacteria. 

Synchronization: synchronization of dynamical systems refers to a process wherein two (or many) systems (either 
identical or nonidcntical) adjust a given property of their motion to a common behavior due to a coupling or to a 
forcing (periodical or noisy). 

Lur'e system: a Lur'e system is a linear dynamic system, feedback interconnected to a static nonlinearity that 
satisfies a sector condition. 



APPENDIX 
A. Notations: 

Throughout this paper, denotes the transpose of a square matrix A. The notation M > (<) is used to define 
a real symmetric positive definite (negative definite) matrix, i?'" denotes the m-dimensional Euclidean space; and 
j^nxm (jgnotes the set of all n x to real matrices. In this paper, if not explicitly stated, matrices are assumed to 
have compatible dimensions. || • || stands for the usual L2 norm of a vector, or the usual spectral norm of a square 
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matrix. Matrix U G Ji^xN jg defined as an irreducible matrices with zeros row sum, whose ofF-diagonal elements are 
all non-positive. The Kronecker product A (g) B of an n x m matrix A and a, p x q matrix B is the np x mq matrix 
defined as 



Ai^B 



AiiB ■ ■ ■ AimB 



AjilB ■ ■ ■ AnmB 



B. Theoretical Analysis of the Synchronization 



Proof of Theorem 1 : We let x{t) = [x{{t), - ■ ■ , xlf{t)\ G i? , and define a Lyapunov function of the following 



form: 



According to 



V{x{t))^x^{t){U ®P)x{t) (11) 
(pp. 136-137), V{x{t)) is equivalent to the following form 

V{x{t)) = ^(-C/y)(x.(t) ~ x,{t)f P{x.{t) - x,{t)). 

i<j 

Hence, if the time derivative of V{x{t)) along the trajectories of (6) is negative, then according to Lyapunov's direct 
method, the genetic oscillators will achieve synchronization. 
Calculating the time derivative of V{x{t)), we have 



V{x{t)) = 2E.<,(-C/.,)(2;.(t) - x,{t)rP[{A - T){xm - + B,{J{xi) - f{x,)) ~ B2{g(x.) - g{x,))] 

+2x'^{t){U ® P)(G ® D + / ® T)xit) 

where T G is an arbitrary real matrix. For all i, j (i ^ j), we have 

= 2{x,{t) - x,{t)yP[{A - T){xdt) - x,{t)) + B,{f{xdt)) - f{x,{t))) - B2{g{x,{t)) - g{x,{t)))] 
< 2{x,{t) - x,{t)fP{A - T){x,{t) - x,{t)) + 2{x,{t) - x,{t)YPBi{f{x,{t)) - f{x,{t))) 
-2{x,{t)-x,{t)fPB2{g{xdt))-g{x,{t))) 

-2Er=i MfiMt)) - fi{x,i{t)))[{fi{xa{t)) - fi{x,i{t))) - ku{xu{t) - x„m 
-'^Y.?=i>'2iigi{xii{t)) - gi{xji{t)))[{gi{xu{t)) - gi{xji{t))) - k2i{xii{t) - x.ji{t))] 
By letting Q = PT, we have 



(12) 
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for alH,j. = [{x,{t) - Xj{t)f , {f{x,{t)) - f{xj{t)))'^ , (dixtit)) - g{xj{t)))'^f G i?3rixi^ Therefore, the first term of 
(12) is negative except for x^{t) = Xj{t), ViJ. By letting Q = PT, {U ® P){G D + 1 (»Tf + {U P){G ® D + 1 (»T) 
is equivalent to the second matrix inequality in (8), which means that the second term of (12) is non-positive. We 
have V{x{t)) < in (12). Thus the Theorem 1 is proved. □ 

Next, we consider the case of coupled nonidentical genetic oscillators (7). In the authors studied the robust 
synchronization of master-slave coupled two nonidentical chaotic systems of the Lur'c form. Here we extend the result 
to the case of complex networks, and apply it to the coupled nonidentical genetic oscillators. 

Proof of Theorem 2: We also use the Lyapunov function (11). By calculating the time derivative of V{x{t)), we 
have 

V{x{t)) - 2E.<,(-f/.,)(a;.(0 + AA,(t)) -T)(x.(t) -a;,(i)) 

+ iBi + ABu{t))if{x,) - f{xj)) - {B2 + AB2.(t))(.9(a;.) - g{x,)) + (AB2« - AB2,)] (13) 
+2x'^{t){U ® P){G ®D + I® T)x{t) 
By letting Q = PT, it is easy to know the second term is nonpositive, and for all «, j (i ^ j), we have 

2{x,{t) ~ Xj{t)YP[{{A + AA,(t)) - T){x,{t) - xj{t)) + (Bi + ABu{mf{x,) - f{x,)) 

-{B2 + AB2^{t)){g{x^) - g{Xj)) + {AB2^ - AB2j)] 

= 2{x,{t) - x,{t) fP[{A - T){x,{t) - x,{t)) + B^{f{x,{t)) - f{x,{t))) - B2{g{x,{t)) - g{x,{t)))] 
+2{x,{t) - x,{t)YP[AA,{t)x,{t) - AA,{t)x,{t) + ABu{t)f{x,{t)) - ABij{t)f{x,{t)) 
~AB2^^t)g{x^{t)) + AB2j{t)g{xj{t)) + AB2^{t) - AB2j{t)] 

< ^,j{t)Mi^,j{t) + 2\\x^{t) ~ X,{t)\\Xrnax{P)i2aidi + 20262 + 20363 + 203) 
= ^^,{t)^U,,{t) - 7||x,(i) - X,{t)f + PWx^it) - X,it)\\ 

< --r\\x,{t) - Xj{t)\\^ + (3\\x,{t) - x,{t)\\ 

where ^.y(t) and Afi are the same as those defined in the Proof of Theorem 1. We have V{x{t)) < ii \\xi{t) — Xj{t)\\ > 
(3/ J. Since 

Xrmn{P)J2i-Uv)\\iXr{t)-Xj{t))f < Vix{t)) < XrnaAP)Y.(-U^Mx^it)~X,{t))f, 
i<j i<j 

we have 

Xm^niP)J2i-U^Mx^it)-X,{t))f < A„,, (P) ^ ^ (-[/,,) . 
i<j i<j 
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Thus we obtain the error bound in Theorem 2. Note that \"""''>p\ is the conditional number of matrix P. □ 
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